In this notebook, we will examine the prior and posterior predictive distributions. This will tell us how well defined our priors are and whether or not they can capture some given data. In addition, we will look at the posterior predictive distribution, to see how well samples from our posterior capture our data after optimisation.
Make sure the environment is correct and print for posterity.
using Pkg
Pkg.status()
Status `~/Projects/NetworkTopology/Project.toml` [b5ca4192] AdvancedVI v0.1.1 [76274a88] Bijectors v0.8.14 [a93c6f00] DataFrames v0.22.5 [2b5f629d] DiffEqBase v6.57.5 `https://github.com/SciML/DiffEqBase.jl.git#sensitivity_interpolation` [41bf760c] DiffEqSensitivity v6.42.0 [0c46a032] DifferentialEquations v6.16.0 [31c24e10] Distributions v0.24.13 [ced4e74d] DistributionsAD v0.6.19 [f6369f11] ForwardDiff v0.10.16 [7073ff75] IJulia v1.23.1 [093fc24a] LightGraphs v1.3.5 [c7f686f2] MCMCChains v4.7.0 [91a5bcdd] Plots v1.10.4 [c3e4b0f8] Pluto v0.12.21 [37e2e3b7] ReverseDiff v1.5.0 [f3b207a7] StatsPlots v0.14.19 [fce5fe82] Turing v0.15.10 [e88e6eb3] Zygote v0.6.3
using Random
using DifferentialEquations
using Turing
using Plots
using StatsPlots
using MCMCChains
using LightGraphs
using Distributions
using Base.Threads
plotly()
Turing.setadbackend(:forwarddiff)
Random.seed!(1);
┌ Info: For saving to png with the Plotly backend PlotlyBase has to be installed. └ @ Plots /home/chaggar/.julia/packages/Plots/6EMd6/src/backends.jl:372
function make_graph(N::Int64, P::Float64)
G = erdos_renyi(N, P)
L = laplacian_matrix(G)
A = adjacency_matrix(G)
return L, A
end
make_graph (generic function with 1 method)
N = 5
P = 0.5
L, A = make_graph(N, P);
The second step of the modelling process will be to define the ODE model. For network diffusion, this is given by:
$$ \frac{d\mathbf{u}}{dt} = -\rho \mathbf{L} \mathbf{u} $$We can set this up as a julia function as follows:
function NetworkDiffusion(u, p, t)
du = -p * L * u
end
NetworkDiffusion (generic function with 1 method)
u0 = rand(5)
p = 2
problem = ODEProblem(NetworkDiffusion, u0, (0.0,1.0), p);
sol = solve(problem, AutoTsit5(Rosenbrock23()), saveat=0.05);
plot(sol)
Now that we have the model, we wish to infer the values from some synthetic data. In this case, we will use points generated with the ODE integrator plus some Gaussian noise.
$$ \mathbf{y} = f(\mathbf{u0}, \rho) + \mathcal{N}(0, \sigma) $$$$ \sigma = 0.02 $$Let's do this and plot the data.
data = Array(sol) + 0.02 * randn(size(Array(sol)))
plot(Array(sol)', ylims=(0,1));
scatter!(data')
Next, we can define a generate model to fit this data. We use the following priors.
$$\sigma \approx \Gamma^{-1}(2, 3)$$$$\rho \approx \mathcal{N}(5,10,[0,10])$$ $$\mathbf{u0} \approx \mathcal{N}(0,2,[0,1])$$
@model function fitode(data, problem)
u_n = size(data)[1]
σ ~ InverseGamma(2, 3) # ~ is the tilde character
ρ ~ truncated(Normal(5,10.0),0.0,10)
u ~ filldist(truncated(Normal(0.5,2.0),0.0,1.0), u_n)
prob = remake(problem, u0=u, p=ρ)
predicted = solve(prob, AutoTsit5(Rosenbrock23()),saveat=0.05)
@threads for i = 1:length(predicted)
data[:,i] ~ MvNormal(predicted[i], σ)
end
return ρ, u
end
fitode (generic function with 1 method)
We can sample directly from the prior distributions and generate data using these values. The results of this will be that can visualise the simulation space of our generative model. It's important for our simulation that we are able to generate samples that can reproduce the data. We should be able to see this from the prior predictive model.
prior_chain = sample(fitode(data,problem), Prior(), 10_000)
chain_array = Array(prior_chain)
plot(Array(sol)[2,:], w=1, legend = false)
for k in 1:500
par = chain_array[rand(1:10_000), 1:7]
resol = solve(remake(problem,u0=par[1:5], p=par[6]),AutoTsit5(Rosenbrock23()),saveat=0.05)
plot!(Array(resol)[2,:], alpha=0.5, color = "#BBBBBB", legend = false)
end
scatter!(data[2,:], legend = false)
Sampling: 100%|█████████████████████████████████████████| Time: 0:00:03
model = fitode(data,problem)
chain = sample(model, NUTS(0.65), MCMCThreads(), 1_000, 10, progress=true);
┌ Info: Found initial step size
│ ϵ = 0.2
└ @ Turing.Inference /home/chaggar/.julia/packages/Turing/XLLTf/src/inference/hmc.jl:188
┌ Info: Found initial step size
│ ϵ = 0.2
└ @ Turing.Inference /home/chaggar/.julia/packages/Turing/XLLTf/src/inference/hmc.jl:188
┌ Info: Found initial step size
│ ϵ = 0.4
└ @ Turing.Inference /home/chaggar/.julia/packages/Turing/XLLTf/src/inference/hmc.jl:188
┌ Info: Found initial step size
│ ϵ = 0.2
└ @ Turing.Inference /home/chaggar/.julia/packages/Turing/XLLTf/src/inference/hmc.jl:188
┌ Info: Found initial step size
│ ϵ = 0.05
└ @ Turing.Inference /home/chaggar/.julia/packages/Turing/XLLTf/src/inference/hmc.jl:188
┌ Info: Found initial step size
│ ϵ = 0.4
└ @ Turing.Inference /home/chaggar/.julia/packages/Turing/XLLTf/src/inference/hmc.jl:188
┌ Info: Found initial step size
│ ϵ = 0.4
└ @ Turing.Inference /home/chaggar/.julia/packages/Turing/XLLTf/src/inference/hmc.jl:188
┌ Info: Found initial step size
│ ϵ = 0.05
└ @ Turing.Inference /home/chaggar/.julia/packages/Turing/XLLTf/src/inference/hmc.jl:188
┌ Info: Found initial step size
│ ϵ = 0.2
└ @ Turing.Inference /home/chaggar/.julia/packages/Turing/XLLTf/src/inference/hmc.jl:188
┌ Info: Found initial step size
│ ϵ = 0.2
└ @ Turing.Inference /home/chaggar/.julia/packages/Turing/XLLTf/src/inference/hmc.jl:188
┌ Warning: The current proposal will be rejected due to numerical error(s).
│ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
└ @ AdvancedHMC /home/chaggar/.julia/packages/AdvancedHMC/MIxdK/src/hamiltonian.jl:47
┌ Warning: The current proposal will be rejected due to numerical error(s).
│ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
└ @ AdvancedHMC /home/chaggar/.julia/packages/AdvancedHMC/MIxdK/src/hamiltonian.jl:47
┌ Warning: The current proposal will be rejected due to numerical error(s).
│ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
└ @ AdvancedHMC /home/chaggar/.julia/packages/AdvancedHMC/MIxdK/src/hamiltonian.jl:47
┌ Warning: The current proposal will be rejected due to numerical error(s).
│ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
└ @ AdvancedHMC /home/chaggar/.julia/packages/AdvancedHMC/MIxdK/src/hamiltonian.jl:47
┌ Warning: The current proposal will be rejected due to numerical error(s).
│ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
└ @ AdvancedHMC /home/chaggar/.julia/packages/AdvancedHMC/MIxdK/src/hamiltonian.jl:47
┌ Warning: The current proposal will be rejected due to numerical error(s).
│ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
└ @ AdvancedHMC /home/chaggar/.julia/packages/AdvancedHMC/MIxdK/src/hamiltonian.jl:47
┌ Warning: The current proposal will be rejected due to numerical error(s).
│ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
└ @ AdvancedHMC /home/chaggar/.julia/packages/AdvancedHMC/MIxdK/src/hamiltonian.jl:47
Sampling (10 threads): 100%|████████████████████████████| Time: 0:00:03
summarize(chain)
Summary Statistics parameters mean std naive_se mcse ess rhat Symbol Float64 Float64 Float64 Float64 Float64 Float64 u[1] 0.1848 0.0092 0.0001 0.0001 6624.8952 1.0003 u[2] 0.0278 0.0187 0.0002 0.0003 4886.4228 1.0000 u[3] 0.1138 0.0250 0.0003 0.0004 4151.8727 1.0023 u[4] 0.5335 0.0195 0.0002 0.0003 4573.5282 1.0007 u[5] 0.4211 0.0196 0.0002 0.0003 4664.3748 1.0011 ρ 2.2972 0.2639 0.0026 0.0037 4776.6836 1.0014 σ 0.0415 0.0036 0.0000 0.0000 6231.2817 1.0007
plot(chain)
Now that we have sampled from the posterior using the NUTS sampler, we can use some of these samples to run simulations and plot the posterior predictive model, similarly as we did for the prior predictive model.
posterior_chain_array = Array(chain)
node = 5
plot(ylims=(0,1))
for k in 1:500
par = posterior_chain_array[rand(1:10_000), 1:7]
resol = solve(remake(problem,u0=par[1:5], p=par[6]),AutoTsit5(Rosenbrock23()),saveat=0.05)
plot!(Array(resol)', alpha=0.1, color = "#BBBBBB", legend = false)
end
scatter!(data', legend = false)
plot!(Array(sol)', w=2, legend = false)
The results show that the posterior distribution has converged well to capture the data for each node time series.
Next, we'll run the same process for a network of size N = 50.
N = 50
P = 0.5
L, A = make_graph(N, P);
u0 = rand(50)
p = 2
problem = ODEProblem(NetworkDiffusion, u0, (0.0,1.0), p);
sol = solve(problem, AutoTsit5(Rosenbrock23()), saveat=0.05);
data = Array(sol) + 0.02 * randn(size(Array(sol)))
plot(Array(sol)', ylims=(0,1));
scatter!(data')
prior_chain = sample(fitode(data,problem), Prior(), 10_000)
Sampling: 100%|█████████████████████████████████████████| Time: 0:00:50
Chains MCMC chain (10000×53×1 Array{Float64,3}):
Iterations = 1:10000
Thinning interval = 1
Chains = 1
Samples per chain = 10000
parameters = u[1], u[2], u[3], u[4], u[5], u[6], u[7], u[8], u[9], u[10], u[11], u[12], u[13], u[14], u[15], u[16], u[17], u[18], u[19], u[20], u[21], u[22], u[23], u[24], u[25], u[26], u[27], u[28], u[29], u[30], u[31], u[32], u[33], u[34], u[35], u[36], u[37], u[38], u[39], u[40], u[41], u[42], u[43], u[44], u[45], u[46], u[47], u[48], u[49], u[50], ρ, σ
internals = lp
Summary Statistics
parameters mean std naive_se mcse ess rhat
Symbol Float64 Float64 Float64 Float64 Float64 Float64
u[1] 0.5034 0.2883 0.0029 0.0031 9333.5639 1.0003
u[2] 0.4961 0.2870 0.0029 0.0031 9419.0500 0.9999
u[3] 0.5020 0.2888 0.0029 0.0028 9923.0826 0.9999
u[4] 0.5001 0.2865 0.0029 0.0029 10176.2693 0.9999
u[5] 0.4994 0.2886 0.0029 0.0028 9858.5310 1.0000
u[6] 0.5047 0.2877 0.0029 0.0029 9938.4932 1.0002
u[7] 0.4971 0.2873 0.0029 0.0032 10037.8616 1.0004
u[8] 0.4997 0.2898 0.0029 0.0033 9847.9247 1.0000
u[9] 0.4961 0.2875 0.0029 0.0027 10002.2412 1.0000
u[10] 0.4988 0.2852 0.0029 0.0029 9788.9466 0.9999
u[11] 0.5016 0.2861 0.0029 0.0030 9882.5158 0.9999
u[12] 0.5027 0.2867 0.0029 0.0028 9515.8846 0.9999
u[13] 0.5020 0.2874 0.0029 0.0027 9449.6027 1.0001
u[14] 0.4991 0.2874 0.0029 0.0033 9645.4124 0.9999
u[15] 0.4974 0.2852 0.0029 0.0027 10565.9311 0.9999
u[16] 0.5008 0.2868 0.0029 0.0025 10665.6138 1.0001
u[17] 0.5005 0.2880 0.0029 0.0027 9798.5500 1.0000
u[18] 0.4990 0.2886 0.0029 0.0030 9944.5766 1.0002
u[19] 0.4990 0.2888 0.0029 0.0031 9694.2273 0.9999
u[20] 0.4994 0.2894 0.0029 0.0029 9990.9275 1.0000
u[21] 0.5016 0.2872 0.0029 0.0026 9692.7716 1.0000
u[22] 0.5019 0.2857 0.0029 0.0028 10027.1702 1.0000
u[23] 0.4986 0.2866 0.0029 0.0025 10104.3037 1.0000
⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮
29 rows omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
u[1] 0.0273 0.2538 0.5026 0.7521 0.9748
u[2] 0.0242 0.2467 0.4973 0.7404 0.9748
u[3] 0.0240 0.2519 0.5051 0.7517 0.9754
u[4] 0.0249 0.2518 0.5050 0.7445 0.9739
u[5] 0.0253 0.2513 0.4980 0.7493 0.9748
u[6] 0.0280 0.2571 0.5015 0.7550 0.9746
u[7] 0.0241 0.2482 0.4924 0.7464 0.9732
u[8] 0.0236 0.2443 0.5041 0.7486 0.9748
u[9] 0.0265 0.2462 0.4944 0.7450 0.9742
u[10] 0.0285 0.2558 0.4971 0.7428 0.9728
u[11] 0.0280 0.2556 0.5041 0.7455 0.9746
u[12] 0.0245 0.2566 0.5030 0.7491 0.9759
u[13] 0.0238 0.2525 0.5036 0.7517 0.9750
u[14] 0.0261 0.2521 0.4995 0.7466 0.9741
u[15] 0.0252 0.2535 0.4935 0.7412 0.9750
u[16] 0.0247 0.2568 0.4986 0.7492 0.9746
u[17] 0.0243 0.2492 0.5003 0.7472 0.9750
u[18] 0.0265 0.2492 0.4969 0.7479 0.9761
u[19] 0.0235 0.2504 0.5023 0.7469 0.9744
u[20] 0.0248 0.2465 0.4999 0.7504 0.9726
u[21] 0.0274 0.2536 0.5025 0.7488 0.9756
u[22] 0.0264 0.2532 0.5047 0.7468 0.9695
u[23] 0.0264 0.2506 0.5021 0.7448 0.9726
⋮ ⋮ ⋮ ⋮ ⋮ ⋮
29 rows omitted
chain_array = Array(prior_chain);
plot(Array(sol)[1,:], w=1, legend = false)
for k in 1:500
par = chain_array[rand(1:10_000), 1:52]
resol = solve(remake(problem,u0=par[1:50], p=par[51]),AutoTsit5(Rosenbrock23()),saveat=0.05)
plot!(Array(resol)[1,:], alpha=0.5, color = "#BBBBBB", legend = false)
end
plot!(Array(sol)[1,:], w=2, legend = false)
scatter!(data[1,:], legend = false)
Plotted above is the time series for a single node (blue), data points (pink) and simulations from the prior (grey). Conversely to our prior predictive model for five nodes, the dynamics here show the model reaches equillibium much faster and almost invariably of starting concentration. Thus, for most data points, i.e. after equillibium, the prior captures the data reasonably well. If one images what the posterior should look like for the equillibium data, it would closely resemble the prior. This lends credence to the idea that the likelihood landscape for this particular data is somewhat flat, since, intuitively, one knows that a flat likelihood means the prior and posterior will be similar, if not the same.